Gaining analytic control of parton showers 



Christian W. Bauer and Frank J. Tackmann 
Ernest Orlando Lawrence Berkeley National Laboratory, University of California, Berkeley, CA 94720 

Parton showers are widely used to generate fully exclusive final states needed to compare theoret- 
ical models to experimental observations. While, in general, parton showers give a good description 
of the experimental data, the precise functional form of the probability distribution underlying the 
event generation is generally not known. The reason is that realistic parton showers are required to 
conserve four-momentum at each vertex. In this paper we investigate in detail how four-momentum 
conservation is enforced in a standard parton shower and why this destroys the analytic control 
of the probability distribution. We show how to modify a parton shower algorithm such that it 
conserves four-momentum at each vertex, but for which the full analytic form of the probability 
distribution is known. We then comment how this analytic control can be used to match matrix 
element calculations with parton showers, and to estimate effects of power corrections and other 
uncertainties in parton showers. 



I. INTRODUCTION 

To analyze experimental data at high energy colliders, 
one needs precise theoretical predictions to compare mea- 
surements against. For such comparisons, Monte Carlo 
event generators that simulate fully exclusive events are 
indispensable tools While it is possible to calculate 
simple observable distributions analytically, in most cases 
a direct comparison of such calculations with the data is 
very difficult, because the experimental analyses have to 
impose a variety of cuts, and detector efficiencies are, in 
general, not uniform over the measured distribution. 

Event generators typically simulate events in three sep- 
arate steps: First, a matrix element generator generates 
an event with few final-state partons, based on full ma- 
trix element calculations, which include all interference 
effects of the quantum field theory. Second, a parton 
shower adds additional partons using classical splitting 
probabilities. Finally, all partons arc turned into hadrons 
according to some QCD model of hadronization. 

In this paper we are concerned with the second step 
in an event generator, the parton shower. In practice it 
is impossible to generate all possible partonic final states 
using matrix element calculations, because the number 
of partons in the final state quickly becomes too large. 
In contrast, parton showers are based on splitting func- 
tions which describe the classical probability that a given 
mother parton splits into two daughter partons. Thus, 
with a certain probability described by the splitting func- 
tions, the parton shower turns a final state containing n 
partons into one with n + 1 partons. The advantage of 
using splitting functions over full matrix elements is that 
this procedure can be iterated to take a simple final state 
with a small number of partons and produce many addi- 
tional partons using a Markov Chain. 

The splitting functions describe the splitting in the 
coUinear limit, where the mother and two daughters 
have large energy and small invariant mass. The par- 
ton shower also resums the leading Sudakov double loga- 
rithms. Recently it was shown that the parton shower 
is formally reproduced as the leading order in an effec- 



tive field theory expansion using soft-coUinear effective 
theory Q, which provides, in principle, a framework for 
systematic improvements. For example, one could at- 
tempt to sum large logarithms beyond the leading order 
or study power corrections. However, almost all theoret- 
ical improvements will necessarily go beyond the level of 
classical splitting probabilities. It is thus very unlikely 
that they can be incorporated into a Markov Chain al- 
gorithm. A typical example are matrix element calcu- 
lations, for which it is hard to directly distribute events 
according to, and which must be explicitly matched with 
parton showers to avoid double counting [1, [^, 0| . 

Other important aspects are theoretical uncertainties 
in the generated distribution ff\. They arise from input 
parameters, like as and quark masses, as well as from 
higher-order power and perturbative corrections. Having 
a reliable estimate of these uncertainties becomes crucial 
at the LHC in searches for New Physics signals where the 
Standard Model background can only be estimated using 
Monte Carlo generators. 

Related to these issues, important practical consider- 
ations must be taken into account as well. The three 
steps described above produce a theoretical distribution 
of events. However, the obtained events still have to be 
run through a detector simulation in order to compare 
them directly with experimental data. It is this addi- 
tional step which consumes by far the most computing 
time in practice. Generating a typical event at the LHC 
before the detector simulation takes only a fraction of a 
second per CPU, whereas propagating an event through 
the detector simulation takes several minutes. In ad- 
dition, it takes of the order of a megabyte to store a 
fully simulated event. The available amount of processing 
power and disk space thus leads to considerable practical 
limitations. For example, it is extremely impractical to 
resimulate the full event set each time the theory makes 
progress, or to simulate separate event sets using dif- 
ferent parameter values in order to estimate theoretical 
uncertainties. 

All of these problems can be solved if the exact distri- 
bution of the generated events is known. If this is the 
case, one can reweight the generated events according to 
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a different theoretical distribution, even after the time- 
consuming detector simulation. This makes it possible 
to reuse existing simulated events, and thus allows one 
with small effort to study theoretical uncertainties and to 
include theoretical improvements whenever they become 
available. 

For this reweighting approach to be possible, one has to 
have control of the precise functional form of the proba- 
bility distribution underlying the event generator. For 
the matrix element generator, this is always the case 
by construction. However, it also requires one to have 
analytic control of the parton shower, which is gener- 
ally not the case for the currently used parton showers. 
The reason is that a realistic parton shower needs to en- 
force four-momentum conservation at each vertex. While 
this is strictly speaking a subleading effect, it typically 
generates cross correlations between different splittings. 
Hence, although the basic probability distribution gov- 
erning a single splitting can be obtained analytically, the 
analytic control over the full distribution is lost in the 
standard parton shower algorithms because of the way 
four-momentum conservation is enforced. 

The main purpose of this paper is to show how the 
analytic control over the full probability distribution can 
be regained, and that in this way the reweighting ap- 
proach becomes feasible. In the next section we review 
the basic parton shower and set up our notation. In 
Sec. mil we take the parton shower of Sherpa 0] , which 
has the same algorithm as Pythia's virtuality-ordered 
parton shower liioi . , as an example to investigate a 
real parton shower algorithm in detail. We then show in 
Sec. lIVI how it can be modified to satisfy four- momentum 
conservation at each vertex, while at the same time re- 
taining the analytic control of the probability distribu- 
tion. In Sec. |V] we discuss how these results can be 
applied in the different contexts described above, and 
Sec. IVII contains our conclusions. 



II. SETUP 
A. Branching Probabilities 

The purpose of this section is to review some of the ba- 
sics of parton showers needed in our discussion and, 
in the process, introduce our notation. To be specific, we 
consider a final-state parton shower with the invariant 
mass as evolution variable. For simplicity, we assume all 
particles to be massless, although particle masses can be 
included in a straightforward way. 

The branching of a mother parton with some energy 
E into two daughter partons is determined by two in- 
dependent kinematic variables, which are chosen to be 
the invariant mass of the mother t (or equivalently the 
total invariant mass of the daughters) and the energy 
splitting z, which determines how the mother's energy 
is distributed between the daughters. Before the mother 
is branched, it is still on shell with t = 0. During the 



branching step the mother is put off shell and given an 
invariant mass t > 0. At the same time the energy split- 
ting z is obtained. 

The single branch probability V{t,z), defined as the 
differential probability for a branching to occur with cer- 
tain values t and z, is given by 

7'(i,z) = /(i,z)exp|- r^\t' ""dz7(i',^')| 

= /(<,z)n(t,W), (1) 

where f (t, z) is the usual Altarelli-Parisi splitting func- 
tion [l3| and H(i,imax) is the well-known Sudakov fac- 
tor |14l | , which resums the leading Sudakov double loga- 
rithms. The exact form of tmax and Zcut in Eq- HI) de- 
pends on the details of the parton shower implementation 
and will be discussed later. 

The algorithm to determine the value of t at which 
the branching occurs is thought of as evolving t from 
some high start value tmax down to some lower value. 
The Sudakov factor H(t,tmax) then corresponds to the 
no-branching probability, i.e., the probability that no 
branching between imax and t occurs. To see this, one 
integrates V{t, z), 

/ dt I dz7'(t,z) = 1-H(ti,t,„ax). (2) 

The left-hand side is the probability for the mother to 
branch somewhere between ti and tmax- Since the total 
probability is unity, n(ii,imax) is the probability that 
the mother does not branch between ti and tmax- 

When the mother is branched, the daughters are still 
on shell with zero invariant mass. In the next step the 
daughters themselves are branched and given non zero 
invariant masses, and after that, their daughters, and so 
on, resulting in a treelike structure as shown in Fig. [1] 
As the invariant mass of each daughter must be less than 
that of its mother, t is always decreasing for consecutive 
branchings. When a branching with t < tcut is obtained, 
the corresponding mother is left unbranched. The parton 
shower terminates once no more branchings with t > tcut 
are found. The value of the cutoff tcut is typically chosen 
to be a low scale of order a few GeV'^. 

Each event produced by the parton shower consists 
of a tree of n branches characterized by the set of vari- 
ables {ti, Zi} = {ti, Zi;t2, Z2', ■ ■ ■ ',tn, Zn\, which can later 
be turned into momentum four-vectors. As discussed in 
the Introduction, we would like to be able to explicitly 
compute the probability for a given event with certain 
values {ti, Zi{ to be generated. Ideally, this probability is 
simply given as the product of individual single branch 
probabilities, schematically, 

PMiU, Zi\) = V{ti,z{) X V{t2, Z2) X • • • . (3) 

If the parton shower satisfies Eq. ^ , the problem of find- 
ing a closed form for Pps({ti, Zi}) reduces to working out 
the exact form of the single branch probability V{t, z). 
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FIG. 1; Diagrammatic representation of a tree of branches. 



the whole branch) hnking the daughters' energies to Eq, 

El = zoEo , Er^{1~ zo)Eo . (4) 

Thus, energy conservation is automatically satisfied, and 
El and Eji are not free variables. 

In the rest frame of the mother, the value of zq is fixed 
in terms oHl and tR, 
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(5) 



Using Eqs. dH) and (O, the magnitude of the daughters' 
three-momenta is 

PcM = El - tL ^ Er - tR ^ -rHio, tL,tRf , (6) 



It is important to point out that, even though the par- 
ton shower itself is only valid in the limit where the in- 
variant mass of each daughter is much smaller than that 
of its mother, we still need to know the exact form of 
^ps({ii, Zi}) for all values of {U, Zi}. In other words, it is 
not sufficient for our purposes to only know Pps({^i, Zi}) 
expanded in the limit where the parton shower is valid. 

Eq. ([3]) would be trivially satisfied if all branchings 
were independent of each other. In general, this is not the 
case for two reasons. First, the branching of each daugh- 
ter depends on the initial conditions set by its mother, 
which implies that Pps({ii, -^i}) depends on the specific 
structure of the tree. However, as long as each branching 
only depends on previous branchings, the total proba- 
bility can still be written as a product of single branch 
probabilities as in Eq. ([3]). 

The second reason Eq. ^ can be violated is more in- 
volved and related to the basic issue of the parton shower 
we wish to address. At the time each branch is generated, 
the corresponding daughters have not yet been branched 
and are still on shell. However, the phase-space limits 
for the branch following from four-momentum conser- 
vation depend on the daughters' final invariant masses. 
Hence, only after both daughters have been branched 
themselves, one can come back and enforce the correct 
phase-space limits on the branch. Usually this step in- 
volves some kinematic reshuffling, which ends up intro- 
ducing a complicated correlation between the daughters' 
branches and thereby violating Eq. ([3]). 



B. Kinematics 

We begin by working out the phase space for a single 
1 — !■ 2 branch. We denote all kinematic quantities (see 
Fig. [T|) related to the mother particle with a subscript 
and those of the left and right daughters with a subscript 
L and R, respectively. We regard the mother's invariant 
mass ^0 and energy Eq as fixed and take the daughters 
to have invariant masses II and tn. The energy splitting 
Zq is regarded as a property of the mother (or rather of 



with the usual phase-space factor 

Hto,tL,tR) = ^^J{to~tL-tRY ^AtLtR. (7) 

to 

To enforce timelike daughters one needs 

tL<El, tR<El, (8) 

which, using Eqs. ([6|) and (O, is equivalent to the usual 
phase-space limit 



tL + VtR.< Vh- 



(9) 



The kinematics in a general frame is obtained by boost- 
ing the above re sults. The magnitude of the boost 
Po — \/l — to /Eq is fixed by Eq and to, while its direction 
can be described by the angle 9o between the boost axis 
and the daughters' three-momenta in the mother's rest 
frame. Since encodes the information how the energies 
of the daughters are boosted relative to the mother's rest 
frame, it effectively determines zq, 



Zq - z™+PqC0s9q 



IPcmI 



1 + 7^-7^ + Pacos9oX{ta,tL,tR) 
to to 



(10) 



and vice versa. 



1 2zo-{l+tL/to-tR/to) 

cos Oo = - w, , , X • (11) 

The phase-space limits in a general frame are a simple 
generalization of the limits in the rest frame. 



tL + VtR < Vta , 



|cos6lo| < 1 . 



(12) 



Using Eq. (fTTj) . the limit on cos 6*0 is equivalent to the zq 
limit 



zo 



2 V to to 



<^X{to,tL,tR), (13) 
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commonly found in parton shower algorithms. 

Finally, we look at a double branch 1 ^ 2 — > 4 where 
each daughter further branches into two on-shell parti- 
cles. In this case, the daughters' energy splittings ^l.a, 
or equivalently 0L,i?j are additional free variables. The 
complete phase space now is just an extension of Eq. (fT2|) . 



VU + Vh^.<VU>, |cos0o|<l, |cos0i,fl,| < 1. (14) 

The limits on z^n equivalent to |cos0L^i^| < 1 are anal- 
ogous to Eq. (fT^ with the daughters' invariant masses 
set to zero. Hence, the complete phase space in terms of 
tL^R and zo,L,R reads 



tL + VtR < Vto , 



U) to 



< 



< 



2 



with 



ft 



(15) 



(16) 



and El^ji as in Eq. (g]). 

Eq. (|15p explicitly shows the problem mentioned at 
the end of the previous section. Initially, zq is generated 
assuming tL,R = 0, but since the limit on zq depends 
on tL and ta, the generated value of zq has to be ad- 
justed after {Il^zl) and (tn^zn) have been determined. 
Changing zq, however, changes El^r and Pl,r- This in 
turn changes the limits on zl and zr, which can render 
their values invalid. In addition, and tji are deter- 
mined independently from one another, so the constraint 
V^L + V^R < V^o can be violated as well. 



III. A STANDARD PARTON SHOWER 

To study a concrete example of a standard parton 
shower, we consider the final-state parton shower of 
Sherpa which employs the same algorithm as the 
Pythia virtuality-ordered parton shower [§, [l^, flit . 
Other algorithms which employ different ordering vari- 
ables can be found in Refs. [H [H, [13, 111 • 



A. The Algorithm 

To be able to enforce four-momentum conservation, 
the parton shower algorithm always branches two sisters 
in pairs. That is, in each iteration it takes an existing 
1^2 branch, consisting of a mother and two unbranched 
daughters, and converts it into a 1 ^ 2 ^ 4 double 
branch by branching both daughters. To do so, the algo- 
rithm proceeds in three steps as depicted in Fig. O 



1. Branch both daughters, each according to ^{t, z). 

2. Shuffle zo^zJ™(zo,tL,ti?j. 

3. Check kinematics in terms of new Zq™: 

(a) If successful, accept daughter branches. 

(b) If failed, evolve daughter with larger t further 
down and return to step 2. 

In step 1, each daughter is branched separately, with val- 
ues for (tL, zl) and (t/f, zji) distributed according to the 
single branch probability T'{t, z). In step 2, zq is changed 
to a new value Zq°^ {zo,tL,tR), which is derived from its 
old value and takes into account the now nonzero values 
of tL,R- In the mother's restframe this shuffling sim- 
ply sets Zq to the correct value Eq. ([5]). In a 
general frame the form of Zq™(zo, ti, tn) is not dictated 
by kinematics anymore, but is usually chosen to satisfy 
Eq. (fT3|) . In step 3, the kinematics are checked, using 
the new value Zq™ . If they are satisfied, the daughter 
branches are accepted. Otherwise, the algorithm takes 
the daughter with the larger t, evolves it further down, 
and goes back to step 2. 

In the remainder of this section we discuss the details 
of this algorithm. We first work out the precise form 
of the single branch probability 7'(t, z) employed by the 
algorithm. We then move to discuss in detail steps 2 and 
3, which implement four- momentum conservation, but as 
one can already see from Fig. [5] and the above discussion, 
introduce a complicated correlation between tL and tn, 
which clearly violates Eq. ([3]). 



B. The Single Branch Probability 

In this section we are only interested in a single 1^2 
branching of a mother into two daughters. This means 
that each of the daughters in the algorithm described 
above acts as the mother now, and similarly, the mother 
and the other daughter in the algorithm now act as 
grandmother and sister, respectively. 

In the first step, the algorithm independently generates 
two sets of values (ti,Zi) and (tn^Zfi) according to the 
single branch probability V{t,z), introduced in Eq. ([1]). 
The precise form of V{t,z) that is actually used in the 
algorithm can be written as 



V{t,z) = /(t,z)n(i,w)e(i max 

1 



t) Oit - teut) 



X e 



2cut(^) — 



(17) 



with the Sudakov factor 
U{ti,t2) = cxpj-^ dt 



t(t) 



-2cut(t) 



dz/(<,z) . (18) 



In Eq. ()17p we explicitly included all kinematic 9 func- 
tions restricting the allowed ranges of t and z. All in- 
formation on the precise form of V(t,z) is encoded in 
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b) Failed, decrease larger t 



FIG. 2: Diagrammatic representation of a standard parton shower algorithm. Solid lines represent off-shell partons with nonzero 
invariant mass, dashed lines unbranched, on-shell partons. 



the integration limits imax Sind Zcut(Oi the splitting 
function f{t,z). Note that V{t,z) is not normalized to 
unity, because it describes the differential distribution in 
{t, z) for an entire 1^2 branch. It does not include the 
probability that the mother parton does not branch, in 
which case z is undefined. However, we can still define 
the differential distribution in t for a single parton, 

v{t) = n(t 

cut 7 ^max 

)5{t)+ JdzV{t,z), (19) 

where the first term is the no-branching probability, and 
V{t) is now properly normalized to unity, 

y'dtT'Ct) =n(teut,tmax) + [l-n(teut,tmax)] = 1- (20) 

The precise form of the splitting function f{t, z) de- 
pends on the specific type of splitting (q — s- qg^ g qq, 
or g^ gg), e.g. 



as{p)CF^ 1 1 
~ 1~ 



27r 



1 



(21) 



The scale at which as is evaluated is, in general, a func- 
tion of t and z. For simplicity we use /j,^ = t/A. Another 
typical choice is /x^ = z(l — z)t. 

We stress that the kinematics relevant to Eq. (fTT]) as- 
sumes that both daughters and the mother's sister are 
on-shell, massless particles. The expressions for imax and 
-Zcut(i) arising from the phase space limits are then 



ZcUt{t) — ^ ; 



(22) 



where i^ini is the initial energy of the mother^, and /3 — 
\J\ — V^ini subsection. 



^ Usually, is given in terms of the grandmother's and Eq. 
One exception is the case when the mother is the final parton 
coming from a hard interaction in the grandmother's rest frame. 
In this case, Eini is chosen to be Eini(t) = v'*o(l + Vo)/2, which 
is the exact result for an on-shell, massless sister. 



In addition to the pure phase space limits, the algo- 
rithm includes several restrictions which modify the form 
of imax and Zcut(i)- First, since the parton showers is or- 
dered in the evolution variable, t is always smaller than 
the initial value tjni where the evolution starts, which is 
usually chosen to be the invariant mass of the grand- 
mother, tini = io- Second, the cutoff on the algorithm, 
i > icut, is realized as a cut on = |ptP > icut/4, 
where pr is the daughters' transverse momentum with 
respect to the mother's flight direction. In terms of t and 
z we have 



t t / 



1n2 
2. 



(23) 



Hence, > tcut/4 implies t > tcut, but with the addi- 
tional restriction 



1 

'-2 



(24) 



Furthermore, parton showers require the opening angle 
between the daughters of subsequent emissions to always 
decrease. This angular ordering ensures that the branch- 
ing is correctly described not only in the coUinear limit, 
where the splitting functions arc derived, but also in the 
soft limit, where the branching is coherent [1, [1^. The 
opening angle •& is given by 



cos "d — 1 — 



2z(l 



(25) 



which allows us to translate the angular ordering cut 
"!? < '^cut, where dent is the opening angle of the previous 
branch, into limits on t and z, 



1 - cos ^?cut 



z — 



< ^ / 1 _ ^ 1 -h COS t^cut 
- 2Y /32i;.2^. l-C0Sl?cut 



(26) 



Putting everything together, the single branch prob- 
ability 7^(t, z) depends on several additional parameters 
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restricting the allowed range of t and z, 

r{t, z) = Vit, z\tinU EinU tcut, ??cut) , (27) 
where the limits on t and z are determined by 

2 



1 — COS 



(28) 



COS l?cut / -, ^cut 

^cut(^) = o nim<^ . 1 - — ,\1 — 

2 [ V f3^Ef 1 - cosi^cut V t 



Note that these limits automatically include the phase- 
space limits, Eq. (j^ . which are reproduced for -dcut — ti" 
and tcut — 0. Also, for t ^ E'^, corresponding to the 
mother's rest frame, Zcut{t) 0, forcing z —f 1/2, as 
required. 



C. Enforcing Momentum Conservation 

In each iteration the algorithm in Fig. [2] starts with a 
value for zq which, as discussed above, was determined 
in a previous iteration assuming t^n = and satisfying 
|2zo ~ 1| ^ 00- In particular, in the mother's rest frame 
it starts with zq = 1/2, whereas after step 1, II and 
tf( have become nonzero and so the correct value is now 
Zq^^ given in Eq. ([5]). Thus, the value of zq has to be 
changed (shuffled), or otherwise the kinematics can never 
be satisfied. 

This shuffling happens in step 2 of the algorithm. 
There are various ways to do this, and Sherpa uses 



Now assume, for instance, that II > tR. Then Zq™ > zq, 
which increases (decreases E^) and decreases Pl (in- 
creases Pr). Using Eq. ([TS]) . it follows that the available 
phases space for z^ shrinks, and hence the value of z^ 
may not be allowed anymore. The same is true for the 
right branch in case > t^- 

For this reason, the algorithm has to explicitly check 
the kinematics, which is done in step 3. In Sherpa this 
is implemented by checking various different kinematical 
constraints arising from four-momentum conservation, all 
of which can be reduced to the constraints in Eq. (|15p . 
If all constraints are satisfied, the daughter branches are 
accepted and the algorithm proceeds with the next itera- 
tion. Otherwise, if any phase-space limit is violated, the 
algorithm picks the daughter with the larger t, generates 
new values (t, z) according to V{t, z) with the previous t 
as tini, and goes back to step 2. 

It should be clear from this discussion that steps 2 
and 3 in the algorithm introduce a complicated cross 
correlation between the probabilities to get certain val- 
ues (iL,Zi) and {tR^ZR), which breaks the factorization 
in Eq. ([3]), and makes it virtually impossible to find a 
closed-form expression for the double branch probability 

P{tL,ZL;tR,ZR). 

To end this section, we note that the shuffling of zq 
is proportional to iL.fl/^o- As the parton shower is for- 
mally only valid for t^ R <^ to, this is formally a power 
suppressed effect. Four- momentum conservation, how- 
ever, is an important power correction that must be taken 
into account in order to obtain realistic events, also be- 
cause in the end the parton shower is used for any values 



z^''"(zo,tL,tR) 
1 

' 2 



^-^ + (2zo-l)A(io,tL,ti?J 

To to 



(29) 



with X{to,tL,tR) given in Eq. 0. For zq = 1/2, Eq. ((29l) 
reduces to Eq. (O, as required. Since the original value of 
Zq satisfies |22o — 1| < Po, the shuffling in Eq. ensures 
that Zq always satisfies the correct phase-space limit, 
Eq. p3)) . for t^ R 0. Nevertheless, four-momentum 
conservation can still be violated in two ways. First, 
it may not be possible at all to find a new physical 
value Zq°™. Namely, t^^R may not satisfy the constraint 



V^L + -s/^i? < \/to in Eq. , which ensures timelike 
daughters and that X{to,tL,tR) in Eq. ([29]) is well de- 
fined. This can happen, because the values of and Ir 
were determined independently from one another, only 
subject to the constraint t^^R < tmax [see Eq. ([28)) ]. 



Second, changing zq 
of the two daughters. 



also changes the energies 



TP , TTTICW new 77' TP , TPnOW / 1 ^HOWNtP 



and accordingly [see Eq. 



now 



l^tL,R/El<^^. 



(30) 



(31) 



IV. THE ANALYTIC PARTON SHOWER 
A. The Analytic Algorithm 

It is the final step in the parton shower, where a gener- 
ated set {tL,ZL) or {IrjZr) can be rejected, which leads 
to the complicated correlation in the double branch prob- 
ability. As explained above, there are two reasons the 
kinematics can fail in step 3. First, it is possible that 
the original values oit^ R violate \/tL + \/tR < y/to, and 
second, the required momentum shuffling in step 2 can 
render the values of z^./j invalid. We now show how both 
of these problems can be dealt with without introducing 
complicated correlations. 

We start by looking at the second problem, for which 
it is instructive to understand in more detail the physical 
picture behind the shuffling in Eq. ([29]) . Looking back at 
Eq. (jlOp . one can think of zq in a general frame as given 
in terms of Iq^l^r and /SoCOsSq. Vice versa, for given 
to,L,R, Zq determines the boost factor /3o cos^o- When zq 
was generated, iL,i? — 0, and Eq. pll)l imphes 

Po cos 6io = 220 - 1 • (32) 

Thus, selecting a value for zq is equivalent to choosing 
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1 . Branch left 
► 



{tL,eL) 




2. Branch right 
► 



{tL,eL) 




3. Ok 



FIG. 3: Diagrammatic representation of the analytic algorithm. Solid lines represent ofT-shell partons with nonzero invariant 
mass, dashed lines unbranched, on-shell partons. The evolution of tL in step 1 starts at tini = to, and that of tn in step 2 at 



the boost factor /3o 0036*0. With this in mind, it is easy 
to understand what the shuffling in Eq. ([^5]) is doing 
physically. It simply holds (3o cos 6*0 fixed at its generated 
value, and recomputes zq using Eq. (jlOp after the daugh- 
ters acquire nonzero invariant masses. In other words, 
the algorithm treats the boost factor Pq cos Oq as the fun- 
damental and Zq as the derived quantity.^ 

In the same way, generating values for zl,r in step 1 
in Fig. [2] can be regarded as generating and fixing values 
for the boost factors 



(3l,B. cos eL,R = '2ZL,R - 1 



(33) 



However, Pl.r are not free quantities, but change with 
Zq as functions of Il^r [see Eqs. ([50)) and (^1]) ]. The 
only way to keep PL,RCOsdLji fixed when zg is shuf- 
fled is to balance the change in Pl,r by a corresponding 
change in cosOl h. In this picture, we now immediately 
see what can go wrong. When shuffling zq, either /3l or 
/3i? might decrease too much, resulting in an unphysical 
value 10036*1 > 1 for either 9l or 6*^. 

This picture also leads to a simple and general solution 
to the problem: Instead of /3cos0, we can just as well 
use the boost angle 9 itself as the fundamental quantity. 
When a value for z is generated, it is translated into a 
value for cos 6, 



cos 9 



2z- 1 

/3 



(34) 



which is held fixed at any later stage in the algorithm. 
The advantage is that the phase-space limits 1003 6*1 < 1 
are completely independent of any other kinematic vari- 
ables and are always satisfied. Holding cos 9 fixed in 
Eq. (|34|) corresponds to an additional shuffle 



'(z,/3,/3-") = i 



1 + (2z- 1 



(35) 



to be applied to zl.r. whenever (3l,r change as a result of 
shuffling Zq. Similarly to Eq. Eq. ([55]) ensures that 
zl,r always satisfy their phase-space constraints. 

Note that Eqs. and ([55]) can be realized as two 
special cases of a generalized shuffle, which follows from 
Eqs. §^ and pUl). 



(36) 



tL 
t 



t 



{2z 



old 



/□new 

^)^Kt,tL,tR) 



where is the value used when was generated, and 
^now jg j^g.^ value. This makes the practical imple- 
mentation into the current algorithms straightforward, 
because all one has to do is to store /3°''* and change the 
shuffling function. 

Since all Eq. ([36)1 does is to hold cos 9 fixed, we can also 
go one step further and directly use (t, 9) to describe the 
individual branches, with the daughters' energies El^r 
always given as functions of Eq, 9q and tg^L ji. In this 
way, any necessity to keep track of zo^l,r and how and 
when they are shuffled is removed, which makes the entire 
algorithm very transparent. 

At this point, the only kinematical check that would be 
left in step 3 is the simple constraint ^ytL + \/tR, < -^/to, 
which cannot be eliminated by a change of variables. 
However, we can recast the correlation it introduces into 
a calculable form by branching the daughters in two sep- 
arate steps, as shown in Fig. O In the first step the left 
daughter is branched, starting the evolution at t-mi = to- 
In the second step the right daughter is branched, start- 
ing the evolution at tini = i\/to ~ \fth)^ ■, which auto- 
matically takes care of the remaining constraint. In this 
way, the double branch probability for the left and right 
branches can be written as the product of two single 
branch probabilities 



^ This statement is true in general, even if zq is shuffled more than 
once, because in the algorithm the value of zq used on the right- 
hand side of Eq. I|29|l is always the originally generated value, 
never one which was already obtained from shuffling. 



P(tL,QL\tR, 9r) = V{tL,9L\tini = h) (37) 

X V[tR,9R\tini = iVui- Vt^f] . 

Since the iteration steps in the generation of an entire 
event are already independent, Eq. ([57)1 allows us to write 
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the total probability to generate a given event as a prod- 
uct of single branch probabilities as in Eq. ([3]), which is 
what we set out to do. 



B. Practical Implementation 

In practice, one has several choices in implementing 
this new analytic algorithm in a real parton shower. To 
eliminate the asymmetry between and tji in Eq. p7p . 
one can randomly choose which daughter of a given 
branch acts as the left daughter and is branched first. Al- 
ternatively, one can always branch the daughter with the 
larger (or smaller) initial energy first. A third choice is to 
generate a test value of t for each daughter, call the larger 
one ii, and then branch the other daughter starting the 
evolution at the smaller oHl and {y/to — ^fth)^ . For this 
choice, the double branch probability still factorizes and 
involves an additional factor of the no-branching proba- 
bility n(to, ii). This last choice may be the most natural 
one from the point of view of a global evolution [20| . 

Once the left daughter has been branched, one also 
has the choice which energy to use as the initial energy 
-Eini for branching the right daughter. One could either 
keep the original energy computed from with = 
or recompute it from Qq with the new value of tz,, where 
again the latter choice seems to be the more natural one. 

We have implemented the changes to the algorithm de- 
scribed above in Sherpa's parton shower. For simplicity, 
we always branch the left daughter first and keep the orig- 
inal energy for when branching the right daughter. 
The only things we had to change then were to separate 
the branching of the two daughters and to change the 
function returning a new value of z to use Eq. p6p . As 
a cross check, we did not remove the kinematic checks 
done in step 3 of the original algorithm and tested that 
with our modifications they indeed never fail (apart from 
extremely rare occasions where the failure is due to nu- 
merical inaccuracies). 



C. Comparison and Numerical Results 

We now discuss the impact the change in the algo- 
rithm has on the generated events. The original algo- 
rithm rejects branches in the third step if the kinematics 
is not satisfied, which leads to a lowering of at least one of 
the invariant masses of the daughters that are branched. 
Hence, compared to the analytic algorithm, where this 
third step is absent, the original algorithm suppresses 
large values of t and enhances low values of t. Since the 
effect on the kinematics from having nonzero t is more 
pronounced for larger i, this relative suppression is ex- 
pected to increase with increasing t. 

To estimate the expected size of this effect, we note 
that shuffling z is a power suppressed effect, of order 
t/io, which one can think of as changing the z limits of 
integration. Thus, upon integration over z, the difference 



in the two algorithms corresponds to a power correction 
to the splitting function f{t) = J dz f{t, z), schematically 



A/(t)=/(t)xo(^) 



(38) 



Furthermore, the single branch probability gener- 
ated by the algorithm always has the form P{t) — 
f{t) exp[J dtf{t)]. The integral of a difference in f{t) 
in the exponent gives rise to an additional finite pertur- 
bative difference. Thus, Eq. ((38|) translates into a change 
in the single branch probability 



AP{t) = P{t) X 



0{-]+Oias) 
to 



(39) 



The appearance of perturbative corrections can also be 
understood in another way. Since the total probability 
P{t) is normalized to unity [see Eq. p^ ]. and power 
corrections must vanish for t — > 0, an increase of P{t) at 
large t via power corrections can only be compensated at 
small t by a decrease of P(t) via perturbative corrections. 

To illustrate the difference between the algorithms at 
the level of a single branching, we generated one million 
double branches starting from to = (91.2 GeV)^ in the 
mother's rest frame, and look at the average of the double 
branch probability 



P{t) — / d^L dzi dtfj dzfl 



6{t - tL) + S{t - tL) 



X P{tL, ZL;tR, zr) . 



(40) 



The results are shown on the left of Fig. |4] for the origi- 
nal Sherpa algorithm [dark (blue) triangles] and the an- 
alytic algorithm [medium (orange) dots]. The solid (or- 
ange) line shows the analytic result for P{t) obtained 
from Eq. ([37)l . As expected, it matches the distribution 
generated with the analytic algorithm. One can also see 
the suppression of the original algorithm at large t. The 
difference between the two algorithms for large t is well 
within the expected size of power corrections Eq. (|39p . 
indicated by the gray shading. 

To show the effect on fully showered events, we gen- 
erated one million events with Sherpa using the original 
and the analytic algorithm. A typical physical observ- 
able is the thrust distribution, da/dT, which measures 
the jettiness of a given event, with T — + 1 for two narrow 
jets and T — > 1/2 for spherical events. On the right of 
Fig. [4] we show the integrated thrust distribution 



dT 



l-T 



da 
dT 



(41) 



obtained with the original algorithm [dark (blue) trian- 
gles] and the analytic algorithm [medium (orange) dots]. 
Here, the gray shading gives an estimate of the size of per- 
turbative corrections. As expected, compared to the orig- 
inal algorithm, the analytic algorithm suppresses small 
values of r, corresponding to branchings at small t, but 
within the estimated size of perturbative corrections. 
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FIG. 4: Comparison of the analytic parton shower [medium (orange) dots] and Sherpa's parton shower [dark (blue) triangles]. 
Left: The function P{t) as defined in Eq. (|40|l . The solid (orange) line shows the result obtained from Eq. (I37|l . and the gray 
shading gives an estimate of the expected size of power corrections. Right: The integrated thrust distribution, where the gray 
shading shows the expected size of finite perturbative corrections. The error bars show the statistical uncertainties. 
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FIG. 5: Pull distribution for the double branch probability 
P(tL,tR) defined in Eq. (|42|l . See text for further explanation. 



As a further check that the analytic algorithm dis- 
tributes according to the known probabihty Eq. (|37p . we 
keep the left and right branches separate and consider 
the double branch probability integrated over zl,r 



P{tL,tR)= dzLdziiP{tL,ZL;tR,ZR) . (42) 



The result of the analytic algorithm agrees well within 
the statistical uncertainties with the analytic expecta- 
tion from Eq. ([57]) . This is illustrated in Fig. O where 
we show the pull distribution for P{tL,tR), i.e. the dif- 
ference between the generated and analytic distributions 
divided by the statistical uncertainty of the generated 
distribution. 



V. REWEIGHTING EVENTS 

We have shown how to construct an analytic par- 
ton shower algorithm which distributes events accord- 
ing to a known probability distribution Pps({ii, Zi}). 
This result opens up the possibility to obtain events 
that are distributed according to some other distribution 
PnewiiU, Zi}) by proceeding in three steps: 

1. Generate events using the analytic algorithm de- 
scribed in this work. 



2. Assign the weight 
each event. 



{{k,z,})/Pps{{U,z,}) to 



3. If desired, unweight the event sample by vetoing 
events according to their relative weights. 

The third step is optional and only needed if one desires 
final events with unit weight. Similarly, the allowable 
size of the weights will depend on the specific applica- 
tion. The power of this approach is twofold. First, it can 
provide a very efficient way to distribute according to 
some distribution Pncw, which may not be possible oth- 
erwise. Second, since the reweighting does not change 
the event kinematics, it can be applied at any later stage 
in the event generation, in particular, after the detector 
simulation. We like to stress again that, for this reweight- 
ing approach to be possible, it is essential to know the 
exact form of Pps, i.e., it would be insufficient to only 
know the leading terms in the power expansion of Pps. 
We now discuss two immediate applications. 



A. Distributing and Matching Matrix Elements 



We first consider the case where P„, 



iviE IS given 



by the differential distribution obtained from full matrix 
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element calculations. Even though it is possible to ob- 
tain the squared amplitude, and thus Pme, for a moder- 
ate number of partons in the final state, it is still quite 
difficult to distribute events according to Pme- The rea- 
son is that Pme has large peaks arising from poles in the 
amplitudes due to on-shell intermediate particles. This 
makes it impossible to employ a simple hit-or-miss algo- 
rithm with random numbers drawn from a fiat distribu- 
tion. Instead one has to rely on numerically inverting the 
integral of Pme over phase space. Since the dimensional- 
ity of phase space increases with each additional parton, 
the required numerical integrations quickly become very 
time-consuming. 

The parton shower describes the same IR physics as 
the full QCD amplitude. It thus contains the same poles 
from on-shell intermediate particles, and the resulting 
weights Pme /Pps in step 2 are expected to be moderate. 
Generating events with a parton shower is several or- 
ders of magnitude faster than the numerical phase-space 
integrations required otherwise. Hence, it should be pos- 
sible to accommodate relative weights even of O(IOO) to 
0(1000) and still achieve a reasonably high efficiency. As 
the events will later be run through a detector simulation, 
step 3 has to be included such that the final events have 
unit weight. One way to think of this is that the parton 
shower acts as a phase-space generator that automati- 
cally contains the correct pole structure of the matrix 
elements. 

From this viewpoint, one could consider a much sim- 
pler version of the Markov Chain algorithm that only 
attempts to capture the underlying pole structure, with- 
out trying to address other important effects (e.g. getting 
the correct soft limit via the angular ordering) that are 
already included in the matrix elements. This would pro- 
vide an alternative and efficient algorithm to distribute 
events according to known matrix element expressions. 

Nevertheless, to obtain realistic, fully exclusive events 
one still has to attach a parton shower to the matrix 
element calculations, which is usually nontrivial due to 
double counting issues. There are a few dedicated algo- 
rithms [1, d, with several implementations [2l|, [l^, l23| 
available to consistently match tree-level matrix elements 
for many partons with parton showers. In addition, sev- 
eral approaches are pursued by now to incorporate matrix 
elements for the first hard emission at next-to-leading or- 
der [13, [2^, [23| . To combine these two separate classes 
of matrix element corrections, some work has been car- 
ried out in Ref. [2S|. In this respect, the SCET-based 
approach of Ref. '2"! seems quite promising. 

Using the analytic parton shower, one can generate 
fully showered events with many partons in the final 
state and then reweight the n hardest emissions to the 
matrix element result for n final-state partons. This 
idea is similar to the older merging method used by 
Pythia to correct the first shower emission 0, [l^, [U, 
In our case, one assigns to each event the weight 
PME{{tk, Zk}) / Pps{{tk, Zk}) instep 2, where k = l,...,n 
numbers the n partons with the largest t. It follows that 



any observable sensitive to the distribution in the n hard- 
est partons and inclusive in all other partons will be de- 
termined by the full matrix element result, while any 
further emission is determined by the parton shower. In 
this way the analytic parton shower not only allows one 
to efficiently distribute matrix elements, but in addition 
provides a simple and powerful tool to match matrix el- 
ements and parton showers. An implementation of this 
result will be given elsewhere. Note that the matching is 
completely determined at the analytic level by the form 
of Pme- In principle, it could be carried out for any n 
and at any order in perturbation theory. 



B. Parton Shower Tuning and Uncertainties 

The reweighting can also be performed after the events 
have been run through a detector simulation, which is 
the most time-consuming part of the event generation. 
Hence, one can obtain sets of events with different un- 
derlying distributions, performing only a single run of 
the detector simulation. Of course, to make maximal 
use of the simulated events, one would now like to have 
0(1) weights and also skip step 3. Therefore, to gener- 
ate events one would still use the best available matrix 
element calculations matched with the analytic parton 
shower, as described above. 

One advantage of using the analytic parton shower is 
that one can update already simulated events at any later 
time to the newest theory or parameters. Furthermore, 
having analytic control over all parameters in the par- 
ton shower allows one to estimate uncertainties arising 
from input parameters like asijnz) or quark masses, and 
from higher-order power and perturbative corrections. In 
all cases Pncw is simply given by Pps computed with a 
different set of parameters. Moreover, one could study 
different scheme choices in the parton shower, provided 
of course one has analytic control over Pps correspond- 
ing to the new scheme. For example, one could study the 
scale at which as is evaluated or even the choice of the 
evolution variable. 

As a specific application, we look at power corrections 
and the tuning of the parton shower. Tuning a par- 
ton shower is crucial to get a good description of the 
experimental data, and one is, in effect, adjusting un- 
known power corrections to fit the data. The analytic 
control over the parton shower gives access to a simple 
and systematic way to tune it to data. We can introduce 
power corrections by adding nonsingular terms to the 
splitting function /(t, z) entering Eq. (jl7p . For example, 
for fq^qg{t,z) as in Eq. ([HI), 



/new (^7 ^) 



27r 



1 1- z 



git) 



(43) 



where g{t) is a nonsingular function of t (in general, g 
could also depend on z). This change will affect the 
branching probability for large t, but will leave branches 
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FIG. 6: The function P{t), as defined in Eq. (gOjl, for the an- 
alytic parton shower [medium (orange) dots], after reweight- 
ing [light (green) diamonds], and for Sherpa's parton shower 
[dashed (blue) line]. Error bars indicate statistical uncertain- 
ties. See text for further explanation. 



with small t nearly unaffected. Hence, adjusting g{t) 
changes the power corrections included in the parton 
shower, which can be used to tune the parton shower. 

To illustrate this, we would like to adjust the ana- 
lytic algorithm such that the average of its double branch 
probability, P{t) [see Eq. (HH])], roughly agrees with the 
original algorithm. We use the simplest possible power 
correction, g{t) — a/t„^ax- As discussed in Sec. llVCi 
the analytic algorithm enhances large t compared to the 
original algorithm, so we need a < 0. In Fig. [6] we 
show the result for P{t) from running the analytic algo- 
rithm [medium (orange) dots] and then reweighting each 
event to a new splitting function with a = —1.5 [light 
(green) diamonds], together with the result from run- 
ning the original Sherpa algorithm [dashed (blue) line]. 
The reweighted result agrees remarkably well with the 
original algorithm, given the simple form of the power 
correction. 

While this example only serves as an illustration, the 
analytic control over the tuning parameters is extremely 
useful. After the detector simulation, tuning the par- 
ton shower to data amounts to parametrizing the power 
corrections g{t) in terms of a few parameters and fitting 
them to data. The uncertainties from power corrections 
can then be estimated by varying the tuning parameters 
in some range. 



to generate final states with a large number of partons. 
They rely on splitting functions, which are derived in 
the collinear or soft limit of the underlying theory. The 
conservation of four-momentum in each step of the al- 
gorithm, although technically a subleading effect, is an 
important ingredient to obtain realistic predictions from 
parton showers to compare to data. 

In standard parton shower algorithms, the implemen- 
tation of four-momentum conservation introduces a cross 
correlation between different branches, such that the fi- 
nal probability to produce a given event is not equal to 
the product of individual branching probabilities. In this 
work we have shown that a few simple modifications of 
an existing algorithm yield an analytic parton shower al- 
gorithm that conserves four-momentum at each vertex, 
but distributes according to a known analytic distribu- 
tion. This makes it possible to generate events with exact 
knowledge of the probability with which each event was 
generated, and to reweight the events to a different dis- 
tribution at any later stage in the event generation. 

We have studied the specific case of a virtuality- 
ordered final-state parton shower. Considering the type 
of modifications, we expect the extension to a corre- 
sponding initial-state parton shower to be straightfor- 
ward. It should also be possible to apply the same ideas 
to parton showers using different ordering variables. 

The analytic parton shower proposed here in conjunc- 
tion with the reweighting approach provides a powerful 
tool for experiment and theory, and in the last section 
we have given two examples of this. First, it facilitates 
the distribution of events according to full matrix ele- 
ments, which is otherwise hindered by the need for time- 
consuming phase-space integrations, and at the same 
time provides a very generic way to match matrix ele- 
ments with the parton shower. Second, generated events 
can be reweighted even after they have been run through 
a detector simulation. This allows one to update simu- 
lated events at any later time. It also greatly simplifies 
the study of higher-order corrections and parameter de- 
pendences in the parton shower as well as the estimation 
of uncertainties arising from these sources. For example, 
it provides a convenient way to tune the parton shower 
to data by simply fitting parameters characterizing power 
corrections to the data. 
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